Description

Libraries

library(here)
library(data.table)
library(magrittr)
library(ggplot2)
library(scales)
library(colorspace)
library(janitor)
theme_set(theme_bw())

1 Making maps with R

Data from the North Carolina data

library(sf)
nameshp <- system.file("shape/nc.shp", package = "sf")
d <- st_read(nameshp, quiet = TRUE)
# transform
d <- st_as_sf(d)

# infant deaths 1974
d$vble <- d$SID74
# infant deaths 1979
d$vble2 <- d$SID79
head(d)
## Simple feature collection with 6 features and 16 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -81.74107 ymin: 36.07282 xmax: -75.77316 ymax: 36.58965
## Geodetic CRS:  NAD27
##    AREA PERIMETER CNTY_ CNTY_ID        NAME  FIPS FIPSNO CRESS_ID BIR74 SID74
## 1 0.114     1.442  1825    1825        Ashe 37009  37009        5  1091     1
## 2 0.061     1.231  1827    1827   Alleghany 37005  37005        3   487     0
## 3 0.143     1.630  1828    1828       Surry 37171  37171       86  3188     5
## 4 0.070     2.968  1831    1831   Currituck 37053  37053       27   508     1
## 5 0.153     2.206  1832    1832 Northampton 37131  37131       66  1421     9
## 6 0.097     1.670  1833    1833    Hertford 37091  37091       46  1452     7
##   NWBIR74 BIR79 SID79 NWBIR79                       geometry vble vble2
## 1      10  1364     0      19 MULTIPOLYGON (((-81.47276 3...    1     0
## 2      10   542     3      12 MULTIPOLYGON (((-81.23989 3...    0     3
## 3     208  3616     6     260 MULTIPOLYGON (((-80.45634 3...    5     6
## 4     123   830     2     145 MULTIPOLYGON (((-76.00897 3...    1     2
## 5    1066  1606     3    1197 MULTIPOLYGON (((-77.21767 3...    9     3
## 6     954  1838     5    1237 MULTIPOLYGON (((-76.74506 3...    7     5

2 ggplot2

library(ggplot2)
library(viridis)
# ggplot(d) +
#   geom_sf(aes(fill = vble)) +
#   scale_fill_viridis() + 
#   theme_bw()
ggplot(d) +
  geom_sf(aes(fill = vble)) +
  scale_fill_continuous_sequential() + 
  theme_bw() +
  theme(
    legend.position = "top"
  )

2.1 plotly

library(plotly)
g <- ggplot(d) +
  geom_sf(aes(fill = vble), color = "gray5", size = 0.1) +
  scale_fill_continuous_sequential() + 
  theme_bw()
ggplotly(g)

3 leaflet

The sf object that we pass to leaflet() needs to have a geographic coordinate reference system (CRS) (EPSG code 4326) indicating latitude and longitude. Here, we use the st_transform() function of sf to transform the data d which has CRS given by EPSG code 4267 to CRS with EPSG code 4326.

st_crs(d)
## Coordinate Reference System:
##   User input: NAD27 
##   wkt:
## GEOGCRS["NAD27",
##     DATUM["North American Datum 1927",
##         ELLIPSOID["Clarke 1866",6378206.4,294.978698213898,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4267]]
library(leaflet)

pal <- colorNumeric(palette = "YlOrRd", domain = d$vble)

l <- leaflet(d) %>% 
  addTiles() %>%
  addPolygons(
    color = "white",
    fillColor = ~ pal(vble),
    fillOpacity = 0.8
  ) %>%
  addLegend(pal = pal,
            values = ~ vble,
            opacity = 0.8)
l
l %>% addMiniMap()

3.1 Save interactive figure as png

# Saves map.html

library(htmlwidgets)
saveWidget(widget = l, file = here::here("results/figures", paste(Sys.Date(), "map.html", sep = "_")))

# Takes a screenshot of the map.html created and saves it as map.png
library(webshot)
# webshot::install_phantomjs()
webshot(url = here::here("results/figures", paste(Sys.Date(), "map.html", sep = "_")), file = here::here("results/figures", paste(Sys.Date(), "map.png", sep = "_")))

4 mapview

library(mapview)
mapview(d, zcol = "vble")
library(RColorBrewer)
pal <- colorRampPalette(brewer.pal(9, "YlOrRd"))
# version 2
mapview(d,
        zcol = "vble",
        map.types = "CartoDB.DarkMatter",
        col.regions = pal)
# version 3
map1 <- mapview(d, zcol = "vble")
leaflet::addMiniMap(map1@map)

5 Side by side plots with mapview

library(leaflet.extras2) # similar to patchwork
m1 <- mapview(d, zcol = "vble")
m2 <- mapview(d, zcol = "vble2")
m1 | m2

6 Synchronized maps with leafsync

m1 <- mapview(d, zcol = "vble")
m2 <- mapview(d, zcol = "vble2")
m <- leafsync::sync(m1, m2)
m
htmltools::save_html(html = m, 
                     file = here::here("results/figures", paste(Sys.Date(), "map-side-by-side.html", sep = "_")))

7 tmap

library(tmap)

# define the mode
tmap_mode("plot") # static
# tmap_mode("view") # interactive

tm_shape(d) + 
  tm_polygons("vble") #+

  # tm_shape(st_centroid(d)) + 
  # tm_dots("vble")

8 Mobility flows with flowmapblue

  • To map mobility data
  • data from cellphone companies
# devtools::install_github("FlowmapBlue/flowmapblue.R")
library(flowmapblue)
# data frame 1
locations <- data.frame(
  id = c(1, 2, 3),
  name = c("New York", "London", "Rio de Janeiro"),
  lat = c(40.713543, 51.507425, -22.906241),
  lon = c(-74.011219, -0.127738, -43.180244)
)
locations
##   id           name       lat        lon
## 1  1       New York  40.71354 -74.011219
## 2  2         London  51.50742  -0.127738
## 3  3 Rio de Janeiro -22.90624 -43.180244
# atadate frame 2
flows <- data.frame(
  origin = c(1, 2, 3, 2, 1, 3),
  dest = c(2, 1, 1, 3, 3 , 2),
  count = c(42, 51, 50, 40, 22, 42)
)
flows
##   origin dest count
## 1      1    2    42
## 2      2    1    51
## 3      3    1    50
## 4      2    3    40
## 5      1    3    22
## 6      3    2    42
# final plot of transitions fromand to the 3 destinations
# flowmapblue(locations = locations, flows = flows, mapboxAccessToken = NULL, clustering = TRUE, darkMode = TRUE, animation = TRUE)
LS0tDQpvdXRwdXQ6IA0KICBib29rZG93bjo6aHRtbF9kb2N1bWVudDI6DQogICAgdGhlbWU6IGZsYXRseQ0KICAgIGNzczogInN0eWxlLmNzcyINCiAgICB0b2M6IHRydWUNCiAgICB0b2NfZGVwdGg6IDMNCiAgICB0b2NfZmxvYXQ6DQogICAgICBjb2xsYXBzZWQ6IHRydWUNCiAgICAgIHNtb290aF9zY3JvbGw6IG5vDQogICAgbnVtYmVyX3NlY3Rpb25zOiB0cnVlDQogICAgY29kZV9kb3dubG9hZDogdHJ1ZQ0KICAgIGNvZGVfZm9sZGluZzogc2hvdw0KdGl0bGU6ICJNYWtpbmcgbWFwcyB3aXRoIFIiDQpkYXRlOiAiU3RhcnQ6IDIyLzAzLzIwMjMiDQphdXRob3I6DQogIG5hbWU6IEpvYW8gTWFsYXRvDQogIGFmZmlsaWF0aW9uOiBJbnN0aXR1dG8gZGUgTWVkaWNpbmEgTW9sZWN1bGFyICYgSW1tdW5lLVN0YXRzDQotLS0NCg0KYGBge3Igc2V0dXAsIGluY2x1ZGU9RkFMU0V9DQprbml0cjo6b3B0c19jaHVuayRzZXQoZWNobyA9IFRSVUUsDQogICAgICAgICAgICAgICAgICAgICAgd2FybmluZyA9IEZBTFNFLA0KICAgICAgICAgICAgICAgICAgICAgIG1lc3NhZ2UgPSBGQUxTRSwgDQogICAgICAgICAgICAgICAgICAgICAgZmlnLmFsaWduID0gImNlbnRlciIsDQogICAgICAgICAgICAgICAgICAgICAgZmlnLmFzcCA9IDAuNjE4LA0KICAgICAgICAgICAgICAgICAgICAgIGZpZy53aWR0aCA9IDEwLA0KICAgICAgICAgICAgICAgICAgICAgIGRwaSA9IDEyMCwgDQogICAgICAgICAgICAgICAgICAgICAgb3V0LndpZHRoID0gIjc1JSIpDQpgYGANCg0KIyBEZXNjcmlwdGlvbiB7LX0NCg0KDQoNCiMgTGlicmFyaWVzIHstfQ0KDQpgYGB7cn0NCmxpYnJhcnkoaGVyZSkNCmxpYnJhcnkoZGF0YS50YWJsZSkNCmxpYnJhcnkobWFncml0dHIpDQpsaWJyYXJ5KGdncGxvdDIpDQpsaWJyYXJ5KHNjYWxlcykNCmxpYnJhcnkoY29sb3JzcGFjZSkNCmxpYnJhcnkoamFuaXRvcikNCnRoZW1lX3NldCh0aGVtZV9idygpKQ0KDQpgYGANCg0KIyBNYWtpbmcgbWFwcyB3aXRoIFINCg0KLSBbbGlua10oaHR0cHM6Ly93d3cucGF1bGFtb3JhZ2EuY29tL2Jvb2stZ2RzLzI1LXNwYXRpYWwtbWFraW5nbWFwcy5odG1sKQ0KDQpEYXRhIGZyb20gdGhlIE5vcnRoIENhcm9saW5hIGRhdGENCg0KYGBge3J9DQpsaWJyYXJ5KHNmKQ0KbmFtZXNocCA8LSBzeXN0ZW0uZmlsZSgic2hhcGUvbmMuc2hwIiwgcGFja2FnZSA9ICJzZiIpDQpkIDwtIHN0X3JlYWQobmFtZXNocCwgcXVpZXQgPSBUUlVFKQ0KIyB0cmFuc2Zvcm0NCmQgPC0gc3RfYXNfc2YoZCkNCg0KIyBpbmZhbnQgZGVhdGhzIDE5NzQNCmQkdmJsZSA8LSBkJFNJRDc0DQojIGluZmFudCBkZWF0aHMgMTk3OQ0KZCR2YmxlMiA8LSBkJFNJRDc5DQpoZWFkKGQpDQpgYGANCg0KIyBnZ3Bsb3QyDQoNCmBgYHtyfQ0KbGlicmFyeShnZ3Bsb3QyKQ0KbGlicmFyeSh2aXJpZGlzKQ0KIyBnZ3Bsb3QoZCkgKw0KIyAgIGdlb21fc2YoYWVzKGZpbGwgPSB2YmxlKSkgKw0KIyAgIHNjYWxlX2ZpbGxfdmlyaWRpcygpICsgDQojICAgdGhlbWVfYncoKQ0KZ2dwbG90KGQpICsNCiAgZ2VvbV9zZihhZXMoZmlsbCA9IHZibGUpKSArDQogIHNjYWxlX2ZpbGxfY29udGludW91c19zZXF1ZW50aWFsKCkgKyANCiAgdGhlbWVfYncoKSArDQogIHRoZW1lKA0KICAgIGxlZ2VuZC5wb3NpdGlvbiA9ICJ0b3AiDQogICkNCmBgYA0KDQojIyBwbG90bHkNCg0KYGBge3J9DQpsaWJyYXJ5KHBsb3RseSkNCmcgPC0gZ2dwbG90KGQpICsNCiAgZ2VvbV9zZihhZXMoZmlsbCA9IHZibGUpLCBjb2xvciA9ICJncmF5NSIsIHNpemUgPSAwLjEpICsNCiAgc2NhbGVfZmlsbF9jb250aW51b3VzX3NlcXVlbnRpYWwoKSArIA0KICB0aGVtZV9idygpDQpnZ3Bsb3RseShnKQ0KYGBgDQoNCg0KIyBsZWFmbGV0DQoNClRoZSBzZiBvYmplY3QgdGhhdCB3ZSBwYXNzIHRvIF9sZWFmbGV0KClfIG5lZWRzIHRvIGhhdmUgYSBnZW9ncmFwaGljIGNvb3JkaW5hdGUgcmVmZXJlbmNlIHN5c3RlbSAoQ1JTKSAoRVBTRyBjb2RlIDQzMjYpIGluZGljYXRpbmcgbGF0aXR1ZGUgYW5kIGxvbmdpdHVkZS4gSGVyZSwgd2UgdXNlIHRoZSBfc3RfdHJhbnNmb3JtKClfIGZ1bmN0aW9uIG9mIGBzZmAgdG8gdHJhbnNmb3JtIHRoZSBkYXRhIGQgd2hpY2ggaGFzIENSUyBnaXZlbiBieSBFUFNHIGNvZGUgNDI2NyB0byBDUlMgd2l0aCBFUFNHIGNvZGUgNDMyNi4NCg0KDQoNCmBgYHtyfQ0Kc3RfY3JzKGQpDQpgYGANCg0KYGBge3J9DQpsaWJyYXJ5KGxlYWZsZXQpDQoNCnBhbCA8LSBjb2xvck51bWVyaWMocGFsZXR0ZSA9ICJZbE9yUmQiLCBkb21haW4gPSBkJHZibGUpDQoNCmwgPC0gbGVhZmxldChkKSAlPiUgDQogIGFkZFRpbGVzKCkgJT4lDQogIGFkZFBvbHlnb25zKA0KICAgIGNvbG9yID0gIndoaXRlIiwNCiAgICBmaWxsQ29sb3IgPSB+IHBhbCh2YmxlKSwNCiAgICBmaWxsT3BhY2l0eSA9IDAuOA0KICApICU+JQ0KICBhZGRMZWdlbmQocGFsID0gcGFsLA0KICAgICAgICAgICAgdmFsdWVzID0gfiB2YmxlLA0KICAgICAgICAgICAgb3BhY2l0eSA9IDAuOCkNCmwNCmBgYA0KDQoNCmBgYHtyfQ0KbCAlPiUgYWRkTWluaU1hcCgpDQpgYGANCg0KIyMgU2F2ZSBpbnRlcmFjdGl2ZSBmaWd1cmUgYXMgcG5nDQoNCmBgYHtyLCBldmFsPUZBTFNFfQ0KIyBTYXZlcyBtYXAuaHRtbA0KDQpsaWJyYXJ5KGh0bWx3aWRnZXRzKQ0Kc2F2ZVdpZGdldCh3aWRnZXQgPSBsLCBmaWxlID0gaGVyZTo6aGVyZSgicmVzdWx0cy9maWd1cmVzIiwgcGFzdGUoU3lzLkRhdGUoKSwgIm1hcC5odG1sIiwgc2VwID0gIl8iKSkpDQoNCiMgVGFrZXMgYSBzY3JlZW5zaG90IG9mIHRoZSBtYXAuaHRtbCBjcmVhdGVkIGFuZCBzYXZlcyBpdCBhcyBtYXAucG5nDQpsaWJyYXJ5KHdlYnNob3QpDQojIHdlYnNob3Q6Omluc3RhbGxfcGhhbnRvbWpzKCkNCndlYnNob3QodXJsID0gaGVyZTo6aGVyZSgicmVzdWx0cy9maWd1cmVzIiwgcGFzdGUoU3lzLkRhdGUoKSwgIm1hcC5odG1sIiwgc2VwID0gIl8iKSksIGZpbGUgPSBoZXJlOjpoZXJlKCJyZXN1bHRzL2ZpZ3VyZXMiLCBwYXN0ZShTeXMuRGF0ZSgpLCAibWFwLnBuZyIsIHNlcCA9ICJfIikpKQ0KYGBgDQoNCg0KIyBtYXB2aWV3DQoNCmBgYHtyfQ0KbGlicmFyeShtYXB2aWV3KQ0KbWFwdmlldyhkLCB6Y29sID0gInZibGUiKQ0KYGBgDQoNCmBgYHtyfQ0KbGlicmFyeShSQ29sb3JCcmV3ZXIpDQpwYWwgPC0gY29sb3JSYW1wUGFsZXR0ZShicmV3ZXIucGFsKDksICJZbE9yUmQiKSkNCiMgdmVyc2lvbiAyDQptYXB2aWV3KGQsDQogICAgICAgIHpjb2wgPSAidmJsZSIsDQogICAgICAgIG1hcC50eXBlcyA9ICJDYXJ0b0RCLkRhcmtNYXR0ZXIiLA0KICAgICAgICBjb2wucmVnaW9ucyA9IHBhbCkNCmBgYA0KDQpgYGB7cn0NCiMgdmVyc2lvbiAzDQptYXAxIDwtIG1hcHZpZXcoZCwgemNvbCA9ICJ2YmxlIikNCmxlYWZsZXQ6OmFkZE1pbmlNYXAobWFwMUBtYXApDQpgYGANCg0KDQoNCiMgU2lkZSBieSBzaWRlIHBsb3RzIHdpdGggbWFwdmlldw0KDQpgYGB7cn0NCmxpYnJhcnkobGVhZmxldC5leHRyYXMyKSAjIHNpbWlsYXIgdG8gcGF0Y2h3b3JrDQptMSA8LSBtYXB2aWV3KGQsIHpjb2wgPSAidmJsZSIpDQptMiA8LSBtYXB2aWV3KGQsIHpjb2wgPSAidmJsZTIiKQ0KbTEgfCBtMg0KYGBgDQoNCg0KIyBTeW5jaHJvbml6ZWQgbWFwcyB3aXRoIGxlYWZzeW5jDQoNCmBgYHtyfQ0KbTEgPC0gbWFwdmlldyhkLCB6Y29sID0gInZibGUiKQ0KbTIgPC0gbWFwdmlldyhkLCB6Y29sID0gInZibGUyIikNCm0gPC0gbGVhZnN5bmM6OnN5bmMobTEsIG0yKQ0KbQ0KYGBgDQoNCg0KYGBge3IsIGV2YWw9RkFMU0V9DQpodG1sdG9vbHM6OnNhdmVfaHRtbChodG1sID0gbSwgDQogICAgICAgICAgICAgICAgICAgICBmaWxlID0gaGVyZTo6aGVyZSgicmVzdWx0cy9maWd1cmVzIiwgcGFzdGUoU3lzLkRhdGUoKSwgIm1hcC1zaWRlLWJ5LXNpZGUuaHRtbCIsIHNlcCA9ICJfIikpKQ0KYGBgDQoNCg0KIyBgdG1hcGANCg0KYGBge3J9DQpsaWJyYXJ5KHRtYXApDQoNCiMgZGVmaW5lIHRoZSBtb2RlDQp0bWFwX21vZGUoInBsb3QiKSAjIHN0YXRpYw0KIyB0bWFwX21vZGUoInZpZXciKSAjIGludGVyYWN0aXZlDQoNCnRtX3NoYXBlKGQpICsgDQogIHRtX3BvbHlnb25zKCJ2YmxlIikgIysNCiAgIyB0bV9zaGFwZShzdF9jZW50cm9pZChkKSkgKyANCiAgIyB0bV9kb3RzKCJ2YmxlIikNCmBgYA0KDQojIE1vYmlsaXR5IGZsb3dzIHdpdGggYGZsb3dtYXBibHVlYA0KDQorIFRvIG1hcCBtb2JpbGl0eSBkYXRhDQorIGRhdGEgZnJvbSBjZWxscGhvbmUgY29tcGFuaWVzDQoNCmBgYHtyfQ0KIyBkZXZ0b29sczo6aW5zdGFsbF9naXRodWIoIkZsb3dtYXBCbHVlL2Zsb3dtYXBibHVlLlIiKQ0KbGlicmFyeShmbG93bWFwYmx1ZSkNCmBgYA0KDQpgYGB7cn0NCiMgZGF0YSBmcmFtZSAxDQpsb2NhdGlvbnMgPC0gZGF0YS5mcmFtZSgNCiAgaWQgPSBjKDEsIDIsIDMpLA0KICBuYW1lID0gYygiTmV3IFlvcmsiLCAiTG9uZG9uIiwgIlJpbyBkZSBKYW5laXJvIiksDQogIGxhdCA9IGMoNDAuNzEzNTQzLCA1MS41MDc0MjUsIC0yMi45MDYyNDEpLA0KICBsb24gPSBjKC03NC4wMTEyMTksIC0wLjEyNzczOCwgLTQzLjE4MDI0NCkNCikNCmxvY2F0aW9ucw0KDQojIGF0YWRhdGUgZnJhbWUgMg0KZmxvd3MgPC0gZGF0YS5mcmFtZSgNCiAgb3JpZ2luID0gYygxLCAyLCAzLCAyLCAxLCAzKSwNCiAgZGVzdCA9IGMoMiwgMSwgMSwgMywgMyAsIDIpLA0KICBjb3VudCA9IGMoNDIsIDUxLCA1MCwgNDAsIDIyLCA0MikNCikNCmZsb3dzDQpgYGANCg0KYGBge3J9DQojIGZpbmFsIHBsb3Qgb2YgdHJhbnNpdGlvbnMgZnJvbWFuZCB0byB0aGUgMyBkZXN0aW5hdGlvbnMNCiMgZmxvd21hcGJsdWUobG9jYXRpb25zID0gbG9jYXRpb25zLCBmbG93cyA9IGZsb3dzLCBtYXBib3hBY2Nlc3NUb2tlbiA9IE5VTEwsIGNsdXN0ZXJpbmcgPSBUUlVFLCBkYXJrTW9kZSA9IFRSVUUsIGFuaW1hdGlvbiA9IFRSVUUpDQpgYGANCg0KDQoNCg==